{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test scipy.signal.correlate on some atl06 data from foundation ice stream"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy, sys, os, pyproj, glob, re, h5py\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.signal import correlate\n",
    "from astropy.time import Time\n",
    "\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test scipy.signal.correlate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate some test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "af544ec368c148678aef606cda1e8657",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'index')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dx = 0.1\n",
    "x = np.arange(0,10,dx)\n",
    "y = np.zeros(np.shape(x))\n",
    "ix0 = 30\n",
    "ix1 = 30 + 15\n",
    "y[ix0:ix1] = 1\n",
    "\n",
    "fig,axs = plt.subplots(1,2)\n",
    "axs[0].plot(x,y,'k')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x)), y,'k')\n",
    "axs[1].set_xlabel('index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a signal to correlate with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "edfdfc8b29a04ad5b4cc1494406f5bac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'black = original, blue = shifted')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imposed_offset = int(14/dx) # 14 meters, in units of samples\n",
    "\n",
    "x_noise = np.arange(0,50,dx) # make the vector we're comparing with much longer\n",
    "y_noise = np.zeros(np.shape(x_noise))\n",
    "y_noise[ix0 + imposed_offset : ix1 + imposed_offset] = 1\n",
    "\n",
    "# uncomment the line below to add noise\n",
    "# y_noise = y_noise * np.random.random(np.shape(y_noise))\n",
    "\n",
    "fig,axs = plt.subplots(1,2)\n",
    "\n",
    "axs[0].plot(x,y,'k')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x)), y, 'k')\n",
    "axs[1].set_xlabel('index')\n",
    "\n",
    "axs[0].plot(x_noise,y_noise, 'b')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x_noise)), y_noise,'b')\n",
    "axs[1].set_xlabel('index')\n",
    "\n",
    "fig.suptitle('black = original, blue = shifted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try scipy.signal.correlate:\n",
    "\n",
    "mode ='full' returns the entire cross correlation; could be 'valid' to return only non- zero-padded part\n",
    "\n",
    "method = direct (not fft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = correlate(y_noise,y, mode = 'full', method = 'direct') \n",
    "norm_val = np.sqrt(np.sum(y_noise**2)*np.sum(y**2))\n",
    "corr = corr / norm_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the dimensions of corr?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "corr:  (599,)\n",
      "x:  (100,)\n",
      "x:  (500,)\n"
     ]
    }
   ],
   "source": [
    "print('corr: ', np.shape(corr))\n",
    "print('x: ', np.shape(x))\n",
    "print('x: ', np.shape(x_noise))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3e99b01d48824f5d92c2e4eae49638b3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'Shift  140  samples, or  14.0  m to line up signals')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# lagvec = np.arange(0,len(x_noise) - len(x) + 1)\n",
    "lagvec = np.arange( -(len(x) - 1), len(x_noise), 1)\n",
    "shift_vec = lagvec * dx\n",
    "\n",
    "ix_peak = np.arange(len(corr))[corr == np.nanmax(corr)][0]\n",
    "best_lag = lagvec[ix_peak]\n",
    "best_shift = shift_vec[ix_peak]\n",
    "\n",
    "fig,axs = plt.subplots(3,1)\n",
    "\n",
    "axs[0].plot(lagvec,corr)\n",
    "axs[0].plot(lagvec[ix_peak],corr[ix_peak], 'r*')\n",
    "axs[0].set_xlabel('lag (samples)')\n",
    "axs[0].set_ylabel('correlation coefficient')\n",
    "\n",
    "axs[1].plot(shift_vec,corr)\n",
    "axs[1].plot(shift_vec[ix_peak],corr[ix_peak], 'r*')\n",
    "axs[1].set_xlabel('shift (m)')\n",
    "axs[1].set_ylabel('correlation coefficient')\n",
    "\n",
    "axs[2].plot(x + best_shift, y,'k')\n",
    "axs[2].plot(x_noise, y_noise, 'b--')\n",
    "axs[2].set_xlabel('shift (m)')\n",
    "\n",
    "\n",
    "fig.suptitle(' '.join(['Shift ', str(best_lag), ' samples, or ', str(best_shift), ' m to line up signals']))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's try with some ATL06 data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load some repeat data:\n",
    "\n",
    "\n",
    "import readers, etc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ! cd ..; [ -d pointCollection ] || git clone https://www.github.com/smithB/pointCollection.git\n",
    "# sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "\n",
    "# !python3 -m pip install --user git+https://github.com/tsutterley/pointCollection.git@pip\n",
    "import pointCollection as pc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "moa_datapath = '/srv/tutorial-data/land_ice_applications/'\n",
    "datapath = '/home/jovyan/shared/surface_velocity/FIS_ATL06/'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example hf5 file, if you need to look at the fields\n",
    "# datapath='/home/jovyan/shared/surface_velocity/FIS_ATL06_small/processed_ATL06_20191129105346_09700511_003_01.h5'\n",
    "# !h5ls -r /home/jovyan/shared/surface_velocity/FIS_ATL06_small/processed_ATL06_20191129105346_09700511_003_01.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geographic setting : Foundation Ice Stream"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/jovyan/.local/lib/python3.7/site-packages/pointCollection/__init__.py\n"
     ]
    }
   ],
   "source": [
    "print(pc.__file__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-76.  -74.5 -74.5 -76.  -76. ]\n",
      "[ -98.  -98. -102. -102.  -98.]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a1808ce5927f42c6bedf83e2bd00ba4a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-1676950., -1495950.,  -352175.,  -213175.]), 'origin': 'lower'}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Mosaic of Antarctica for Pine Island Glacier')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# something wrong with pointCollection\n",
    "\n",
    "spatial_extent = np.array([-102, -76, -98, -74.5])\n",
    "lat=spatial_extent[[1, 3, 3, 1, 1]]\n",
    "lon=spatial_extent[[2, 2, 0, 0, 2]]\n",
    "print(lat)\n",
    "print(lon)\n",
    "# project the coordinates to Antarctic polar stereographic\n",
    "xy=np.array(pyproj.Proj(3031)(lon, lat))\n",
    "# get the bounds of the projected coordinates \n",
    "XR=[np.nanmin(xy[0,:]), np.nanmax(xy[0,:])]\n",
    "YR=[np.nanmin(xy[1,:]), np.nanmax(xy[1,:])]\n",
    "MOA=pc.grid.data().from_geotif(os.path.join(moa_datapath, 'MOA','moa_2009_1km.tif'), bounds=[XR, YR])\n",
    "\n",
    "# show the mosaic:\n",
    "plt.figure()\n",
    "MOA.show(cmap='gray', clim=[14000, 17000])\n",
    "plt.plot(xy[0,:], xy[1,:])\n",
    "plt.title('Mosaic of Antarctica for Pine Island Glacier')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load repeat track data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ATL06 reader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    # below: file_re = regular expression, it will pull apart the regular expression to get the information from the filename\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "        D['data_start_utc'] = h5f['/ancillary_data/data_start_utc'][:]\n",
    "        D['delta_time'] = h5f['/' + beam + '/land_ice_segments/delta_time'][:]\n",
    "        D['segment_id'] = h5f['/' + beam + '/land_ice_segments/segment_id'][:]\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in files; this next cell took ~1 minute early in the morning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read 10 data files of which 0 gave errors\n"
     ]
    }
   ],
   "source": [
    "# find all the files in the directory:\n",
    "# ATL06_files=glob.glob(os.path.join(datapath, 'PIG_ATL06', '*.h5'))\n",
    "\n",
    "ATL06_files=glob.glob(os.path.join(datapath, '*.h5'))\n",
    "\n",
    "D_dict={}\n",
    "error_count=0\n",
    "for file in ATL06_files[:10]:\n",
    "    try:\n",
    "        D_dict[file]=atl06_to_dict(file, '/gt2l', index=slice(0, -1, 25), epsg=3031)\n",
    "    except KeyError as e:\n",
    "        print(f'file {file} encountered error {e}')\n",
    "        error_count += 1\n",
    "print(f\"read {len(D_dict)} data files of which {error_count} gave errors\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot ground tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "390d3abc578849b4aa9c4bb5cefca894",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-1676950., -1495950.,  -352175.,  -213175.]), 'origin': 'lower'}\n",
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-1676950., -1495950.,  -352175.,  -213175.]), 'origin': 'lower'}\n"
     ]
    }
   ],
   "source": [
    "plt.figure(figsize=[8,8])\n",
    "hax0=plt.gcf().add_subplot(211, aspect='equal')\n",
    "MOA.show(ax=hax0, cmap='gray', clim=[14000, 17000]);\n",
    "hax1=plt.gcf().add_subplot(212, aspect='equal', sharex=hax0, sharey=hax0)\n",
    "MOA.show(ax=hax1, cmap='gray', clim=[14000, 17000]);\n",
    "for fname, Di in D_dict.items():\n",
    "    cycle=Di['cycle']\n",
    "    if cycle <= 2:\n",
    "        ax=hax0\n",
    "    else:\n",
    "        ax=hax1\n",
    "    #print(fname)\n",
    "    #print(f'\\t{rgt}, {cycle}, {region}')\n",
    "    ax.plot(Di['x'], Di['y'])\n",
    "    if True:\n",
    "        try:\n",
    "            if cycle  < 3:\n",
    "                ax.text(Di['x'][0], Di['y'][0], f\"rgt={Di['rgt']}, cyc={cycle}\", clip_on=True)\n",
    "            elif cycle==3:\n",
    "                ax.text(Di['x'][0], Di['y'][0], f\"rgt={Di['rgt']}, cyc={cycle}+\", clip_on=True)\n",
    "        except IndexError:\n",
    "            pass\n",
    "hax0.set_title('cycles 1 and 2');\n",
    "hax1.set_title('cycle 3+');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Map view elevations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ad662fad0196446cb79b8ab818088b59",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "map_fig=plt.figure()\n",
    "map_ax=map_fig.add_subplot(111)\n",
    "# MOA.show(ax=map_ax, cmap='gray', clim=[14000, 17000])\n",
    "for fname, Di in D_dict.items():\n",
    "    # select elevations with good quality_summary\n",
    "    good=Di['atl06_quality_summary']==0\n",
    "    ms=map_ax.scatter( Di['x'][good], Di['y'][good],  2, c=Di['h_li'][good], \\\n",
    "                  vmin=0, vmax=1000, label=fname)\n",
    "map_ax._aspect='equal'\n",
    "plt.colorbar(ms, label='elevation');\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat track elevation profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ben Smiths's code to plot the individual segments:\n",
    "def plot_segs(D6, ind=None, **kwargs):\n",
    "    \"\"\"\n",
    "    Plot a sloping line for each ATL06 segment\n",
    "    \"\"\"\n",
    "    if ind is None:\n",
    "        ind=np.ones_like(D6['h_li'], dtype=bool)\n",
    "    #define the heights of the segment endpoints.  Leave a row of NaNs so that the endpoints don't get joined\n",
    "    h_ep=np.zeros([3, D6['h_li'][ind].size])+np.NaN\n",
    "    h_ep[0, :]=D6['h_li'][ind]-D6['dh_fit_dx'][ind]*20\n",
    "    h_ep[1, :]=D6['h_li'][ind]+D6['dh_fit_dx'][ind]*20\n",
    "    # define the x coordinates of the segment endpoints\n",
    "    x_ep=np.zeros([3,D6['h_li'][ind].size])+np.NaN\n",
    "    x_ep[0, :]=D6['x_atc'][ind]-20\n",
    "    x_ep[1, :]=D6['x_atc'][ind]+20\n",
    "\n",
    "    plt.plot(x_ep.T.ravel(), h_ep.T.ravel(), **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A revised code to plot the elevations of segment midpoints (h_li):\n",
    "def plot_elevation(D6, ind=None, **kwargs):\n",
    "    \"\"\"\n",
    "    Plot midpoint elevation for each ATL06 segment\n",
    "    \"\"\"\n",
    "    if ind is None:\n",
    "        ind=np.ones_like(D6['h_li'], dtype=bool)\n",
    "    # pull out heights of segment midpoints\n",
    "    h_li = D6['h_li'][ind]\n",
    "    # pull out along track x coordinates of segment midpoints\n",
    "    x_atc = D6['x_atc'][ind]\n",
    "\n",
    "    plt.plot(x_atc, h_li, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7615572623514a019756611bb6179f02",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D_2l={}\n",
    "D_2r={}\n",
    "\n",
    "# specify the rgt here:\n",
    "rgt=\"0027\"\n",
    "rgt=\"0848\" #Ben's suggestion\n",
    "\n",
    "# iterate over the repeat cycles\n",
    "for cycle in ['03','04','05','06','07']:\n",
    "    for filename in glob.glob(os.path.join(datapath, f'*ATL06_*_{rgt}{cycle}*_003*.h5')):\n",
    "        try:\n",
    "            # read the left-beam data\n",
    "            D_2l[filename]=atl06_to_dict(filename,'/gt2l', index=None, epsg=3031)\n",
    "            # read the right-beam data\n",
    "            D_2r[filename]=atl06_to_dict(filename,'/gt2r', index=None, epsg=3031)\n",
    "            # plot the locations in the previous plot\n",
    "            map_ax.plot(D_2r[filename]['x'], D_2r[filename]['y'],'k');  \n",
    "            map_ax.plot(D_2l[filename]['x'], D_2l[filename]['y'],'k');\n",
    "        except Exception as e:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "plt.figure();\n",
    "for filename, Di in D_2l.items():\n",
    "    #Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "    hl=plot_elevation(Di, ind=Di['atl06_quality_summary']==0, label=f\"cycle={Di['cycle']}\")\n",
    "    #hl=plt.plot(Di['x_atc'][Di['atl06_quality_summary']==0], Di['h_li'][Di['atl06_quality_summary']==0], '.', label=f\"cycle={Di['cycle']}\")\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pull out a segment and cross correlate: \n",
    "\n",
    "Let's try x_atc = 2.935e7 thru 2.93e7 (just from looking through data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "08a85087275443f79286e3b26e34ab51",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "cycles = [] # names of cycles with data\n",
    "for filename, Di in D_2l.items():\n",
    "    cycles += [str(Di['cycle']).zfill(2)]\n",
    "cycles.sort()\n",
    "    \n",
    "# x1 = 2.93e7\n",
    "# x2 = 2.935e7\n",
    "\n",
    "beams = ['gt1l','gt1r','gt2l','gt2r','gt3l','gt3r']\n",
    "\n",
    "### extract and plot data from all available cycles\n",
    "fig, axs = plt.subplots(4,1)\n",
    "x_atc = {}\n",
    "h_li = {}\n",
    "h_li_diff = {}\n",
    "times = {}\n",
    "for cycle in cycles:\n",
    "    # find Di that matches cycle:\n",
    "    Di = {}\n",
    "    x_atc[cycle] = {}\n",
    "    h_li[cycle] = {}\n",
    "    h_li_diff[cycle] = {}\n",
    "    times[cycle] = {}\n",
    "\n",
    "    filenames = glob.glob(os.path.join(datapath, f'*ATL06_*_{rgt}{cycle}*_003*.h5'))\n",
    "    for filename in filenames:\n",
    "        try:\n",
    "            for beam in beams:\n",
    "                Di[filename]=atl06_to_dict(filename,'/'+ beam, index=None, epsg=3031)\n",
    "\n",
    "                times[cycle][beam] = Di[filename]['data_start_utc']\n",
    "                \n",
    "                # extract h_li and x_atc for that section\n",
    "                x_atc_tmp = Di[filename]['x_atc']\n",
    "                h_li_tmp = Di[filename]['h_li']#[ixs]\n",
    "                \n",
    "                # segment ids:\n",
    "                seg_ids = Di[filename]['segment_id']\n",
    "#                 print(len(seg_ids), len(x_atc_tmp))\n",
    "                \n",
    "                # make a monotonically increasing x vector\n",
    "                # assumes dx = 20 exactly, so be carefull referencing back\n",
    "                ind = seg_ids - np.nanmin(seg_ids) # indices starting at zero, using the segment_id field, so any skipped segment will be kept in correct location\n",
    "                x_full = np.arange(np.max(ind)+1) * 20 + x_atc_tmp[0]\n",
    "                h_full = np.zeros(np.max(ind)+1) + np.NaN\n",
    "                h_full[ind] = h_li_tmp\n",
    "                \n",
    "                \n",
    "                x_atc[cycle][beam] = x_full\n",
    "                h_li[cycle][beam] = h_full\n",
    "                                   \n",
    "                             \n",
    "#                 ### here is where you would put a filter\n",
    "#                 # you would want to de-mean and detrend that section first:\n",
    "#                 h = h_full\n",
    "#                 x = x_full\n",
    "#                 h = h - np.nanmean(h) # de-mean\n",
    "#                 h = scipy.signal.detrend(h, type = 'linear') # de-trend; need to deal with nans first\n",
    "#                 # use scipy.signal.filter to filter\n",
    "\n",
    "#                 # differentiate that section of data\n",
    "                h_diff = (h_full[1:] - h_full[0:-1]) / (x_full[1:] - x_full[0:-1])\n",
    "                h_li_diff[cycle][beam] = h_diff\n",
    "\n",
    "                # plot\n",
    "                axs[0].plot(x_full, h_full)\n",
    "                axs[1].plot(x_full[1:], h_diff)\n",
    "#                 axs[2].plot(x_atc_tmp[1:] - x_atc_tmp[:-1])\n",
    "                axs[2].plot(np.isnan(h_full))\n",
    "                axs[3].plot(seg_ids[1:]- seg_ids[:-1])\n",
    "\n",
    "\n",
    "\n",
    "        except:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0f3cc08a85a941479e6dcf85238e2b3b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e7c32a0965c54b3ab5ee79955e5836ea",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "921b81a0da124a5d9212fdfe354bc3db",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7d0cd544cb814743968299fa4ecedde6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd018861e2394da09ea02dbf8f07a0ff",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a64ebdc0ef024133961bbd51b7a89276",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_veloc = len(cycles) - 1\n",
    "\n",
    "segment_length = 5000 # m\n",
    "x1 = 2.925e7#x_atc[cycles[0]][beams[0]][1000] <-- the very first x value in a file; doesn't work, I think b/c nans # 2.93e7\n",
    "\n",
    "search_width = 1000 # m\n",
    "dx = 20 # meters between x_atc points\n",
    "\n",
    "for veloc_number in range(n_veloc):\n",
    "    cycle1 = cycles[veloc_number]\n",
    "    cycle2 = cycles[veloc_number+1]\n",
    "    t1_string = times[cycle1]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "    t1 = Time(t1_string)\n",
    "                       \n",
    "    t2_string = times[cycle2]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "    t2 = Time(t2_string)\n",
    "    \n",
    "    dt = (t2 - t1).jd # difference in julian days\n",
    "        \n",
    "    velocities = {}     \n",
    "    for beam in beams:\n",
    "        fig1, axs = plt.subplots(4,1)\n",
    "        \n",
    "        # cut out small chunk of data at time t1 (first cycle)\n",
    "        x_full_t1 = x_atc[cycle1][beam]\n",
    "        ix_x1 = np.arange(len(x_full_t1))[x_full_t1 >= x1][0]\n",
    "        ix_x2 = ix_x1 + int(np.round(segment_length/dx))      \n",
    "        x_t1 = x_full_t1[ix_x1:ix_x2]\n",
    "        h_li1 = h_li_diff[cycle1][beam][ix_x1-1:ix_x2-1] # start 1 index earlier because the data are differentiated\n",
    "        \n",
    "        # cut out a wider chunk of data at time t2 (second cycle)\n",
    "        x_full_t2 = x_atc[cycle2][beam]\n",
    "        ix_x3 = ix_x1 - int(np.round(search_width/dx)) # offset on earlier end by # indices in search_width\n",
    "        ix_x4 = ix_x2 + int(np.round(search_width/dx)) # offset on later end by # indices in search_width\n",
    "        x_t2 = x_full_t2[ix_x3:ix_x4]\n",
    "        h_li2 = h_li_diff[cycle2][beam][ix_x3:ix_x4]\n",
    "\n",
    "        # plot data\n",
    "        axs[0].plot(x_t2, h_li2, 'r')\n",
    "        axs[0].plot(x_t1, h_li1, 'k')\n",
    "        axs[0].set_xlabel('x_atc (m)')\n",
    "        \n",
    "        # correlate old with newer data\n",
    "        corr = correlate(h_li1, h_li2, mode = 'valid', method = 'direct') \n",
    "        norm_val = np.sqrt(np.sum(h_li1**2)*np.sum(h_li2**2)) # normalize so values range between 0 and 1\n",
    "        corr = corr / norm_val\n",
    "        \n",
    "        \n",
    "#         lagvec = np.arange( -(len(h_li1) - 1), len(h_li2), 1)# for mode = 'full'\n",
    "#         lagvec = np.arange( -int(search_width/dx) - 1, int(search_width/dx) +1, 1) # for mode = 'valid'\n",
    "        lagvec = np.arange(- int(np.round(search_width/dx)), int(search_width/dx) +1,1)# for mode = 'valid'\n",
    "\n",
    "        shift_vec = lagvec * dx\n",
    "\n",
    "        ix_peak = np.arange(len(corr))[corr == np.nanmax(corr)][0]\n",
    "        best_lag = lagvec[ix_peak]\n",
    "        best_shift = shift_vec[ix_peak]\n",
    "        velocities[beam] = best_shift/(dt/365)\n",
    "\n",
    "        axs[1].plot(lagvec,corr)\n",
    "        axs[1].plot(lagvec[ix_peak],corr[ix_peak], 'r*')\n",
    "        axs[1].set_xlabel('lag (samples)')\n",
    "\n",
    "        axs[2].plot(shift_vec,corr)\n",
    "        axs[2].plot(shift_vec[ix_peak],corr[ix_peak], 'r*')\n",
    "        axs[2].set_xlabel('shift (m)')\n",
    "\n",
    "        # plot shifted data\n",
    "        axs[3].plot(x_t2, h_li2, 'r')\n",
    "        axs[3].plot(x_t1 - best_shift, h_li1, 'k')\n",
    "        axs[3].set_xlabel('x_atc (m)')\n",
    "        \n",
    "        axs[0].text(x_t2[100], 0.6*np.nanmax(h_li2), beam)\n",
    "        axs[1].text(lagvec[5], 0.6*np.nanmax(corr), 'best lag: ' + str(best_lag) + '; corr val: ' + str(np.round(corr[ix_peak],3)))\n",
    "        axs[2].text(shift_vec[5], 0.6*np.nanmax(corr), 'best shift: ' + str(best_shift) + ' m'+ '; corr val: ' + str(np.round(corr[ix_peak],3)))\n",
    "        axs[2].text(shift_vec[5], 0.3*np.nanmax(corr), 'veloc of ' + str(np.round(best_shift/(dt/365),1)) + ' m/yr')\n",
    "\n",
    "        \n",
    "    fig1.suptitle('black = older cycle data, red = newer cycle data to search across')\n",
    "\n",
    "\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_veloc = len(cycles) - 1\n",
    "\n",
    "segment_length = 1000 # m\n",
    "\n",
    "search_width = 500 # m\n",
    "dx = 20 # meters between x_atc points\n",
    "\n",
    "x1 = 2.915e7#x_atc[cycles[0]][beams[0]][1000] <-- the very first x value in a file; doesn't work, I think b/c nans # 2.93e7\n",
    "x1s = x_atc[cycles[veloc_number]][beams[0]][search_width:-segment_length-2*search_width:10]\n",
    "velocities = {}     \n",
    "for beam in beams:\n",
    "    velocities[beam] = np.empty_like(x1s)\n",
    "for xi,x1 in enumerate(x1s):\n",
    "    for veloc_number in range(n_veloc):\n",
    "        cycle1 = cycles[veloc_number]\n",
    "        cycle2 = cycles[veloc_number+1]\n",
    "        t1_string = times[cycle1]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "        t1 = Time(t1_string)\n",
    "\n",
    "        t2_string = times[cycle2]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "        t2 = Time(t2_string)\n",
    "\n",
    "        dt = (t2 - t1).jd # difference in julian days\n",
    "\n",
    "        for beam in beams:\n",
    "            # cut out small chunk of data at time t1 (first cycle)\n",
    "            x_full_t1 = x_atc[cycle1][beam]\n",
    "            ix_x1 = np.arange(len(x_full_t1))[x_full_t1 >= x1][0]\n",
    "            ix_x2 = ix_x1 + int(np.round(segment_length/dx))      \n",
    "            x_t1 = x_full_t1[ix_x1:ix_x2]\n",
    "            h_li1 = h_li_diff[cycle1][beam][ix_x1-1:ix_x2-1] # start 1 index earlier because the data are differentiated\n",
    "\n",
    "            # cut out a wider chunk of data at time t2 (second cycle)\n",
    "            x_full_t2 = x_atc[cycle2][beam]\n",
    "            ix_x3 = ix_x1 - int(np.round(search_width/dx)) # offset on earlier end by # indices in search_width\n",
    "            ix_x4 = ix_x2 + int(np.round(search_width/dx)) # offset on later end by # indices in search_width\n",
    "            x_t2 = x_full_t2[ix_x3:ix_x4]\n",
    "            h_li2 = h_li_diff[cycle2][beam][ix_x3:ix_x4]\n",
    "\n",
    "            # correlate old with newer data\n",
    "            corr = correlate(h_li1, h_li2, mode = 'valid', method = 'direct') \n",
    "            norm_val = np.sqrt(np.sum(h_li1**2)*np.sum(h_li2**2)) # normalize so values range between 0 and 1\n",
    "            corr = corr / norm_val\n",
    "\n",
    "\n",
    "    #         lagvec = np.arange( -(len(h_li1) - 1), len(h_li2), 1)# for mode = 'full'\n",
    "    #         lagvec = np.arange( -int(search_width/dx) - 1, int(search_width/dx) +1, 1) # for mode = 'valid'\n",
    "            lagvec = np.arange(- int(np.round(search_width/dx)), int(search_width/dx) +1,1)# for mode = 'valid'\n",
    "\n",
    "            shift_vec = lagvec * dx\n",
    "            \n",
    "            if all(np.isnan(corr)):\n",
    "                velocities[beam][xi] = np.nan\n",
    "            else:\n",
    "                ix_peak = np.arange(len(corr))[corr == np.nanmax(corr)][0]\n",
    "                best_lag = lagvec[ix_peak]\n",
    "                best_shift = shift_vec[ix_peak]\n",
    "                velocities[beam][xi] = best_shift/(dt/365)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f79ee35d5e94430688af18e1e218b05d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f0d2fc99910>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "ax1 = plt.subplot(211)\n",
    "for filename, Di in D_2l.items():\n",
    "    #Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "    hl=plot_elevation(Di, ind=Di['atl06_quality_summary']==0, label=f\"cycle={Di['cycle']}\")\n",
    "    #hl=plt.plot(Di['x_atc'][Di['atl06_quality_summary']==0], Di['h_li'][Di['atl06_quality_summary']==0], '.', label=f\"cycle={Di['cycle']}\")\n",
    "plt.legend()\n",
    "plt.ylabel('elevation');\n",
    "\n",
    "ax2 = plt.subplot(212,sharex=ax1)\n",
    "for beam in beams:\n",
    "    plt.plot(x1s+dx*(segment_length/2),velocities[beam],'.',alpha=0.2,ms=3,label=beam)\n",
    "plt.ylabel('velocity (m/yr)')\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylim(0,1500)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:notebook] *",
   "language": "python",
   "name": "conda-env-notebook-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
